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We present a mean field model for coagulation {A + A A) and annihilation {A + A 0) 
reactions on lattices of traps with a distribution of depths reflected in a distribution of mean escape 
times. The escape time from each trap is exponentially distributed about the mean for that trap, 
and the distribution of mean escape times is a power law. Even in the absence of reactions, the 
distribution of particles over sites changes with time as particles are caught in ever deeper traps, 
that is, the distribution exhibits aging. Our main goal is to explore whether the reactions lead to 
further (time dependent) changes in this distribution. 

PACS numbers: 02.50.Ey,82.40.-g,82.33.-z,05.90.+m 

I. INTRODUCTION 

Chemical reactions limited by the motions of the reactants are abundant in nature, and among the most thoroughly 
studied are diffusion-limited reactions. A ubiquitous approach to these systems simply adds the diffusion and reaction 
contributions together so as to construct appropriate reaction-diffusion equations. It is implicit and even explicit in 
these approaches that the diffusive component describes a motion without a memory, and if one invokes an underlying 
continuous time random walk (CTRW) where walkers react when they meet, it is understood that the waiting time 
distributions for reactants to remain at one location before moving on have a finite mean. The most frequently 
invoked waiting time in this scenario is exponential. It is also understood that in fact many microscopic models can 
be subsumed under the same mesoscopic reaction-diffusion umbrella [J 0) B S Q • 

On the other hand, contrary to the diffusive case, it is by now fairly clear that different microscopic scenarios of 
reactions among subdiffusive species, even simple scenarios, lead to different mesoscopic descriptions [gI ItI. IsI. [qI. ITo| . 
For instance, a popular mesoscopic vehicle, the CTRW, involves waiting time distributions of ensembles of particles 
undergoing reactions. The forms of these distributions depend on the underlying microscopic rules and may vary 
from one microscopic scenario to another. It is thus risky to simply assume a form for these distributions as one 
would for diffusive particles. Instead, a more detailed derivation starting from a set of microscopic rules to arrive at 
a mesoscopic CTRW description is necessary. The specific microscopic reaction scenario of interest to us is a lattice 
whose sites are occupied by traps of varying depths, that is, the quenched trap scenario. If independent particles 
simply walk on this lattice, their escape time from each trap is exponentially distributed, with the distribution of 
trap depths reflected in a distribution of mean escape times. An average over this distribution leads to a CTRW 
model in which all sites have the same waiting time distribution . This latter process is spatially homogeneous 

and semi-Markovian, and anomalous diffusion arises if the resulting mean waiting time on each site diverges. This 
is the so-called annealed trap scenario. However, if particles can also react with one another, interesting questions 
immediately arise. Of interest to us are the reactions A -I- A — > A (coagulation) and A -I- A (annihilation); some 
uncertainties arising from different microscopic descriptions can be illustrated with the coagulation reaction. Suppose 
there is an A at a site, and a second A arrives. In the underlying spatially disordered quenched trap model, it does 
not matter which particle is "killed" by the reaction since the waiting time distribution for particle departure from 
any site is exponential. Such a process is often called memoryless because the future of each particle is determined 
only by its present state and not its past, and therefore the differences between the two victims become irrelevant. 
However if, as we subsequently do, we implement a mean field assumption which turns out to lead to a description 
that is asymptotically equivalent to a CTRW with diverging mean waiting times, that is, to an annealed trap model, 
it may make a kinetic difference which of the two particles is "killed" in each reaction. It is not clear in general how 
to include such fine distinctions in a mesoscopic formalism, nor is it clear what sort of rule at the mesoscopic or mean 
field level best mimics the behavior of the underlying trapping problem. 

Our mean field model shares with other CTRW models with diverging mean waiting times the phenomenon of 
"aging" [H, [M, 15, 16]. Even in non-reactive systems, aging causes the waiting time distribution itself to change with 
time as particles settle into sites with ever longer waiting times (corresponding to particles caught in ever deeper traps 
in the underlying system) . Aging for non-reactive particles does not occur in a given finite quenched trap environment 
since the distribution of particles approaches a Boltzmann distribution ■ The most intriguing unanswered question 
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that we address in this work is whether the reactions cause a change in the waiting time distribution. We anticipate 
our answer to this latter question: we find that in our mean field model the reaction does not cause an additional 
change in the waiting time distribution. 

We will compare our model predictions with the results of numerical simulations of the actual distributed trap 
scenario. We find that the mean field approach works well for higher dimensions {d > 3) but not for low dimensions 
(d = 1), a result that agrees with the well-known differences in subdiffusion exponents predicted by CTRW theories 
and those obtained by numerical simulations in the absence of reactions fl§\ The concentrations of surviving particles 
are known to be well reproduced in all dimensions by invoking relations with the number of distinct sites visited in 
the asymptotically equivalent CTRW [1, [23|- On the other hand, the mean field formalism, while restricted to higher 
dimensions, provides additional direct insights into the more detailed information contained in the time dependence 
of the waiting time distributions. 

In Sec. |TT] we describe our mean field model and arrive at a master equation for the density of particles with given 
mean rate for leaving any site in the absence of reactions. This equation explicitly shows the effects of aging. In 
Sec, mil we establish the corresponding master equations in the presence of the coagulation and annihilation reactions. 
In Sec. llVl we discuss the solution of the master equation for the time-dependent rate distribution. Comparisons with 
numerical simulation results for the underlying random trap model are presented in Sec. [V] We conclude with a brief 
summary in Sec. IVII 

II. MEAN-FIELD APPROACH TO A TRAP MODEL 

To construct a mean field model, we start with a lattice whose sites are traps of varying depths. Our model is 
equivalent to that of but our notation is suitably modified to facilitate the inclusion of reactions, which they do 
not consider. The waiting time for leaving a trap i is exponentially distributed, pit\T) = exp(— t/ri), where r^, 
the mean sojourn time in the trap, is determined by the trap's depth via the Kramers (Arrhenius) law. We further 
assume that the distribution of mean waiting times is a power law, for example, p{t) ~ ^t~^~^9{t — 1). Note that 
the asymptotics of the pdf i^{t) of waiting times t upon averaging over this distribution of mean waiting times is then 
of power law form even though the waiting time for each trap is exponentially distributed, i.e., 

/>oo 

^(t)=/ p(t|T)p(T)dT-7r(l+7)T-l-T. (1) 

Jo 

The power law distribution of mean waiting times Tj leads to a power law distribution of the "leaving rates" 
uji ~ \/Ti for departing from a site, 

= 7wT'-i6l(l - w). (2) 

Note that all the moments of the distribution of the rates are finite even if those of the waiting time distribution are 
not. 

Particles are distributed over sites with different mean waiting times, and as time proceeds, the distribution of 
particles over these sites changes even in the absence of reactions because more and more particles get stuck in deeper 
and deeper traps. Our central question concerns the effects of reactions on this evolving distribution. In particular, 
we ask whether the distribution is modified by the reactions. This is difficult to answer, at least analytically, without 
further approximation. The system with distributed traps of varying depths is spatially inhomogeneous, and if the 
particle executes the usual nearest neighbor random walk, the future of any particle moving over this landscape may 
be strongly dependent on its particular location. Furthermore, there may be strong correlations among subsequent 
steps (especially in lower dimensions) because the particle may revisit a previously visited site. To make progress we 
implement a mean field approximation designed to provide information about the evolving distribution of particles over 
the inhomogeneous landscape. Specifically, we assume that the particles do not perform a nearest neighbor random 
walk but instead that each particle is equally likely (probability to step on any of the N sites of the system. In 

other words, the lattice is a complete graph. Equivalently, as an alternative way of viewing the model, we can think 
of particles performing nearest neighbor random walks, but before each step the mean waiting time associated with 
the trap to which the particle is about to step is chosen anew from the distribution p{t) or, correspondingly, a new 
leaving rate is chosen from p{lo), independently of the waiting times chosen in prior steps. It is also equivalent to a 
nearest neighbor CTRW in a very high-dimensional quenched medium. In any case, in the absence of reactions this 
results in a space-homogeneous CTRW model with waiting time distribution ipit) as given in Eq. ([T]) chosen anew at 
each step. 

The assumptions underlying the mean field approach are expected to be more adequate for random walks that are 
transient, that is, ones in which already visited sites are revisited only with a small probability so that almost all 
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sites reached by the walker are new sites. This is the case for random walks in dimensions d> 3 and, as we show in 
Sec. |Vl we do find excellent agreement between the mean field approach and numerical simulations for the reactions 
in quenched trap environments for d > 3. If the walk is recurrent, that is, if the same sites are revisited repeatedly, 
the approximation may be poor. Onc-dimensional walks are recurrent, and d = 2 is the marginal dimension for this 
property. 

As noted earlier, as time evolves the distribution of particles over the sites with given mean waiting times (leaving 
rates) changes even in the absence of the reactions because more and more particles get stuck in deeper and deeper 
traps. We discuss these changes first in the absence and then in the presence of the reactions. For this purpose we 
note that under the mean field assumption, the master equation for the probability ai(t) to be at site i at time t reads 



W = IJ^'^j'^ji*) ~ ^^a^{t). (3) 



d 

di''"^' N 



Since each site is characterized by its own oj, we can pass from the occupation probability to the probability a{Lj,t) 
and introduce the density n{Lu,t) of particles occupying sites with leaving rate between uj and to + duj. One easily 
deduces that n{u},t)duj = a{u!,t)p{u!)duj, where N'p{Lj)duj is the number of sites with leaving rates between ui and 
LD + duj. Note that in a large system 

— (jjjaj(t) ~ — j ujNp{uj)a{uj,t)duj — J ujTi{uj,t)duj. (4) 
The last integral can be associated with the time-dependent mean rate. 



n{t) = J ujn{uj,t)duj. (5) 



Multiplying both sides of Eq. ^ by p{lo) leads to the master equation for n{uj,t), 

—n{uj,t) — ^{t)p{uj) — ujn{uj^t). (6) 

This equation governs the change in the distribution of particles over jump rates in the absence of reactions. Note 
that the total concentration of particles c — J n{u!, t)doj is constant in time. 

III. RATE DISTRIBUTION WITH REACTIONS 

Next we explore the effects of reactions on the site occupation probabilities. We assume that the usual law of mass 
action is appropriate at the local level. For the case of the A + A ^ A reaction within our mean-field approximation, 
the change in the occupation probability at site i is 



d_ 



since the number of particles at a site already occupied by a particle [with probability ai{t)] does not change upon 
the arrival of the new particle. In the case of the A -I- A ^ reaction the number of particles at an occupied site is 
reduced by one upon the arrival of a new particle, so that the corresponding equation reads 

^a,(t) = [l-2a,{t)]^J2^jajit)-cjMt)- (8) 

Again we can focus instead on n{uj,t), but this is no longer a proper probability density because the total number of 
particles is not conserved. Regrouping terms one obtains the reaction equations 



—n{u),t)=p{uj) I UJn(u),t)dio ■ 







ui + LLjn{u!, t)duo 







for the A -f A ^ A reaction, and 



d 
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dt 



UJ + 2 ujn{uj, t)du! 



n{LU,t) (9) 



n{LO,t) (10) 
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for A + A— >0. At t = the particles are homogeneously distributed over all sites in the system, so the initial 
condition for n is n{u;,0) = p{uj). A normalized probability density is obtained by noting that the overall time- 
dependent reactant concentration is given by 

/•OO 

c{t) = / n{uj,t)duj. (11) 

The properly normalized probability density of particles occupying sites with leaving rate between uj and to + duj is 
then given by 

n(uj,t) 

Pi-^t)^^. (12) 
We again introduce the time-dependent mean jumping rate, which is now given by 

/•OO poo 

Q{t)—c^^{t) / ujn{uj,t)duj — / ujp{uj,t)dui, (13) 



and rewrite Eqs. ([9]) and ([T0|) in the form 

d 



^^n(w, t) = c{t)n{t)p{Lu) -[uj + c{t)n{t)]n{LO, t) (14) 
for the A + A ^ A reaction, and 

^n(w, t) = c{t)n{t)p{Lj) - + 2c{t)fl{t)]n{Lj, t) (15) 
for the A + A ^ reaction. Integrating these two equations over the cj-domain gives the classical kinetic equations 

j^c{t) ^ -pmc'it) (16) 

with the stoichiometric coefhcient ( "molarity" ) of the reaction — 1 ioi the A -|- A — > A reaction and /i = 2 for 
A + A ^ 0. The mean jump rate il{t) is thus the time-dependent reaction rate. 

The equation for the probability density p{uj,t) — n{uj,t)c~^{t) corresponding to Eqs. ([H]) and (fT5|) is then 

^v{u,t) +p{oj,t)^^c{t) = n{t)p{uj)-Lop{oj,t) ~ ^icit)n{t)piu;,t). (17) 
dt c[t) dt 

With Eq. we see that the term with the time derivative dc/dt on the left side of this equation and the last term 
on the right side cancel. Thus, the final equation for p{uj, t) does not depend on the molarity /i of the reaction and is 
the same as the equation for n{uj,t) in the absence of reactions, 

j^p{LO,t)^nit)p{Lu)-LOp{u,t). (18) 

This statement also leads to the remarkable conclusion that within the mean field approximation adopted in this 
model the reaction does not affect the waiting time distribution. This answers one of our main questions. 

IV. SOLUTION FOR THE TIME-DEPENDENT RATE DISTRIBUTION 

The principal question posed earlier, namely, whether the reaction changes the waiting time distribution, has been 
answered in the negative within our mean-field approach. It is now instructive to obtain an explicit expression for 
the time dependent reaction rate. Equation is an integro-differential equation for p{u!,t) since fl{t) as defined in 
Eq. (jl3p depends onp{uj^t) itself. However, contrary to the equations for n(aj,t), this equation is linear and can best 
be approached via Laplace transforms. First we consider the time-Laplace transforms of p(uj,t) and ^l{t), 



p(w,s)= / p{uj,t)e-''*dt, n{s) = / n{t)e-'^dt. 
Jo Jo 



(19) 



5 



The transform of Eq. ([18)) is 

sp{uj, s) — p{lu) = fl{s)p{Lu) — ujp{uj, s), (20) 

where we have exphcitly used the initial condition p{uj,0) — p{uj). The formal "solution" for p(aj,s) (with p still 
contained in Q) is 

p{^,s)^[i + n{s)]P^. (21) 

Transforming the definition (fT5|) of n{t) we obtain 

/•OO 

n{s) = / p{uj,s)ujduj, (22) 



so that multiplying both sides of Eq. (|2ip by lu and integrating, we arrive at a closed algebraic equation for fl{s), 



The solution of this equation is 



where 



n{s) ^ [1 + n{s)] f ^^^^dcu. (23) 

Jo s + ^ 



n{s) = ^^^} , (24) 
^ ' l-/(s) 



/(,) ^ r ^d.. (25) 

The integral representing I{s) can be evaluated asymptotically for any long-tailed distribution p(t) using the 
asymptotic method for integrals with weak singularity (see Ch. 1 A§4 in ^l|), as follows. We can rewrite I{s) as 

^i-sT P^cL: ^l-sf P^dc - s r P^du (26) 
Jo s + u; J„ s + uj s + uj 

for any e. Relevant to the long-time asymptotic behavior is the small-s behavior of I{s). As long as e is chosen such 
that e/s— >cxDass— >0, the third integral on the right hand side is of 0(1) in this limit. Furthermore, if p{uj) 'yuj'^~^ 
as LU ~^ and e <C 1, we can write 

I{s) - 1 - S7 / dcj - 0{s) with £ < 1. (27) 

s^o Jq s + uj 

Finally, a change of variables z = e/s then allows us to write for s ~* 

I[s)^l-^s'< dz-0{s) 

Jo ^ + z (28) 

~ 1 ; S '. 

sm 7r7 

For the particular form p{uj) — — uj) introduced earlier, one finds the result valid for all s [Abramowitz and 

Stegun formula 15.3.1] 

i{s)= 2 ^ 2^1(1,1 + 7,2 + 7,-1/^), (29) 
(1 + 7).s 

which for small s (relevant to long-time asymptotic behavior) leads to the second line of Eq. The inverse 

transform of fi(s) then follows upon application of the Tauberian theorem, 

„ , , sin 7r7 „ , , ^ 
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We have thus provided a mesoscopic theoretical foundation for the widely accepted result that the reaction rate decays 
with time. 

The time dependent reaction rate can now be inserted in the rate equation (|16p to solve for the concentration as a 
function of time. At long times we find the explicit result 

fi sm TTj 

We can in fact calculate the full rate distribution by integrating Eq. (HH]). The result up to quadrature is 



p{uj,t) =p{uj)e' 



(32) 



From this, we can extract the u dependence for short times cot ^ 1 as p{u!,t) ~ and for long times ojt 3> 1 as 
p{oj,t) -^p{u})/uj. 

Finally, we note the interesting connections between these results and the number S{t) of distinct sites visited up 
to time t by a particle in a CTRW [l9|. We can connect this quantity to the average reaction rate V,{t) for d > 3 
as follows. Since at each step of the process a particle mostly visits a new site in the system (the random walk is 
"transient" or "non-recurrent" ) , the reaction rate can be approximated by the time derivative of the number of newly 
visited sites, il{t) = dS/dt, and we have [l^ 

si- ip{s) 

where R is the probability of return to the origin (i? = 0.3405 . . . for a simple cubic lattice). Note that the small-s 
behavior of ?A(s) corresponding to the asymptotic behavior of ip{t) in Eq. ^ is 

^(s) ^ 1 _ -J^s^ = 1 _ (s/xy (34) 
sm nj 

Inserting this expression into ((33|) . inserting ([28]) into (I24p . and comparing the resulting expressions, one sees that 
S{s) ^ (1 — R)ri{s)/s for small s, so that ^l{t) ~ (1 — R)dS/dt for large t. Although the interpretation of the reaction 
rate in terms of the number of distinct sites visited is quite standard, the fact that the broad distribution of trapping 
times does not introduce any additional fluctuation effects into the kinetics is not at all trivial. 

A further connection with the distinct number of sites visited in a CTRW occurs for the concentration of surviving 
reactant, namely, c{t) ^ 1/S'(t), a connection that holds not only for d > 2 but also for d < 2, i.e., even when the 
random walk is recurrent 3, 20j. In [l^ we explicitly obtained the 3c? result S{t) ~ [(1 — i?)/r(l + j)]{Xt)'' . In one 
dimension S{t) ^ [\/2/r(l + ^/2)]{Xt)''/^ . We will test the proposition that c{t) ^ i/S{t) in the next section along 
with results of the theory we have developed above. 



V. NUMERICAL RESULTS 



Our discussion supports the expectation that the results of the mean field model approximate those of the underlying 
random depth trap system in Sd (but not in Id). Even in 3d, where walks are transient, it is nevertheless the case that 
there is a finite probability R of return to a previously visited site, e.g. about 1/3 for a simple cubic lattice. In this 
section we compare our results with those of numerical simulations of 3d and Id lattices with traps of random mean 
exit times as described earlier. The algorithm used in the simulations is described in some detail in the Appendix. 

Figure [T] shows a collection of results for the concentration of reactants as a function of time in 3d. The results fall 
essentially on two straight lines, except for finite size effects. The upper set of results is for 7 — 0.5 and the lower 
set for 7 — 0.8. The concentrations shown are indicated as ci{t), which we associate with the coagulation reaction 
A + A ^ A, and C2{t) for the annihilation reaction A + A — > 0. We plot ci{t) and 2c2{t) to ascertain that the 
surviving concentration in the coagulation reaction is twice that of the annihilation reaction. The coincidence of the 
simulation results for the two cases shows that this is indeed the case. Simulation results are shown for two sizes of 
the L X L X L simple cubic lattice with L = 30 and L = 60. The results for the two sizes are the same except for 
the larger value of 7, where we see some deviations from the straight line at very long times due to finite size effects. 
These deviations are more pronounced for the larger 7 (in the faster walk the ends of the lattices are reached earlier) 
and for the smaller lattice. Since this is a log-log plot the straight-line behavior confirms the power law decay of 
the concentrations. Furthermore, we find that the slopes of these lines are close to 7. In particular, we find slopes 
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0.491(2) for 7 = 0.5 and 0.763(2) for 7 = 0.8 (in both cases for both reactfons). The agreement in the former case is 
very good, but somewhat less so in the latter, where we are not in the fully asymptotic scenario [2^ . 

Figure [2] shows our simulation results in Id, again for 7 — 0.5 and 7 = 0.8. The lines in this case are power 
law fits to the simulation results. The size of the lattice is L = 10, 000. The mean field approach is not valid here 
as an approximation to the underlying trap model, and yet a number of interesting results are nevertheless worth 
mentioning. First, we note that the relation ci{t) — 2c2{t) still holds. We do not know whether this indicates that 
here again the reactions do not affect the dynamics of the moving species, but this result indicates that it may be so. 
Second, we note that the exponents of the concentration, while quite different from those of the mean field model, 
agree with those associated with the distinct number of sites visited for the original trap model, namely, that in one 

dimension c(t) ~ i/S{t) ~ (x^) ^ ^-7/(1+7) because all sites within the span of the random walk are visited at 
least once. The numerical fits for the simulation results yield ci{t) — 2c2(t) ~ 0.881(1)^"'^ '^^^®*^^^ for 7 = 0.5 and 
ci{t) = 2c2(i) ~ 0.959(4)^^" '*^°'^' for 7 — 0.8. The exponents 0.33 and 0.43 are in good agreement with the values of 
7/(1 + 7) = 0.333 and 0.444 respectively. 

Finally, we compare our analytic mean field predictions in more detail with numerical simulations of the quenched 
trap model. In Fig. [3] we show results for the distribution p{uj,t) vs lo. The solid lines in both panels correspond to 
the mean field result ([5^ for 7 = 0.5. From high to low on the right side of each panel the curves are for t = 10^, 10^, 
10"*, and 10^ Equation §^ leads to p{uj,t) ~ p{u>) ~ w^-^ for < 1 and p{iJ,t) ^ p{lo)/uj ~ w^-^ for ujt » 1; the 
slopes —0.5 and —1.5 in the log-log plot are evident. The symbols show the corresponding simulation results. The 
top panel shows the results for a 3d quenched trap lattice with nearest neighbor steps, and the bottom panel those 
of of a quenched trap lattice in which steps to any site are equally likely. Our only adjustment in these figures is the 
normalization, which can not be obtained properly from the simulations since they do not cover an infinite range of 
u)S. We have made this adjustment by scaling the simulation results so that the quantity R{t) = /^'".°°' P{uj,t)duj < 1 
evaluated using the scaled results agrees with the one obtained by means of the theoretical result ((52| . Here CLimin 
and Umax define the range of ujs covered by the simulations. These results confirm that already in 3d the mean field 
theory captures the behavior as well as does the actual numerical realization of the model. We have also ascertained 
the agreement between simulations and mean field theory for other values of the subdiffusive exponent 7. 



10"' 
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FIG. 1: The reactant concentrations Ci(t) and 2c2{i) in three dimensions as functions of time for 7 = 0.5 and 7 = 0.8 and 
two lattice sizes. The symbols denote numerical simulation results. The solid lines result from the mean field theory, and the 
dashed lines from the connection with the number of distinct sites visited. 



VI. CONCLUSIONS 

We have presented a mean field theory of for coagulation (A + A ^ A) and annihilation (A + A — > 0) reactions 
on a lattice whose sites are occupied by traps of varying depths. The escape times from these traps are distributed 
exponentially about a mean time, and the distribution of mean escape times is of power law form. We calculate 
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FIG. 2: The reactant concentrations ci(t) and 2c2(t) in one dimension as functions of time for 7 = 0.5 and 7 = 0.8. The 
symbols denote the results of numerical simulations, and the lines are linear fits to these results. 
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0) 

FIG. 3: For all results in this figure, 7 = 0.5. The solid lines are calculated from the analytic result Eq. (|32[1 and the symbols 
are the simulation results for lattices with quenched traps scaled as described in the text. From high to low on the right side of 
each panel the curves are for t = 10^, 10^, 10'', and 10^. Top panel: Sd lattice; bottom panel: lattice with equally likely jumps 
to any site, that is, a complete graph. 

reactant concentrations as a function of time as well as the time-dependent distribution of particles over sites with 
different escape rates. The mean field model is designed with the particular goal of studying this evolving distribution 
and the effects of the reactions on it, and is expected to do well in dimensions d > 3 but not in Id. The noteworthy 
outcome of the model is that this distribution is not changed by the occurrence of the reactions. One consequence of 
this result is that the concentration of surviving reactant in the A + A ^ A reaction is double that of the A-\- A ^ Q 
reactions, that is, ci{t) = 2c2(i). While the model does not shed light on the Id case, numerical simulations here also 
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show that ci{t) = 2c2{t) and thus a reaction- insensitive aging distribution of particles is not ruled out. Numerical 
simulations in d = 3 with quenched traps agree with these predictions. 

While the waiting time distribution is unaffected by the reaction in this mean field model, it is interesting to 
speculate about the assumptions that would have to be made in a spatially translationally invariant model with long- 
tailed waiting time distributions to arrive at this conclusion in a formulation that includes an explicit description of the 
reaction (we remind the reader that this is not an issue in the underlying random trap depth model with exponentially 
distributed waiting times) . The result would seem to be automatically correct for the A + A ^ reaction since every 
reacting pair involves two particles that have not been previously involved in a reaction. However, for the A + A ^ A 
reaction the situation is different. Here one of the reaction partners continues its walk and may participate in a 
later reaction. It may make a difference whether the survivor is the A that arrived at the site immediately before the 
reaction ( "no kill" scenario) , or the A that was there already ( "kill" scenario) , or a choice of one or the other according 
to some probability. While it might be tempting to assume that a random choice of one or the other (e.g. with equal 
probability) would lead to the mean field result obtained above, it is not immediately evident that this choice provides 
exactly the correct compensatory effect. Interestingly, in our numerical simulations we find no statistically significant 
difference between results obtained in the two scenarios. 
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APPENDIX: DESCRIPTION OF THE NUMERICAL ALGORITHM 

In this appendix we describe the steps used in our numerical algorithm. 

1. We generate a Id or 3(i lattice. 

2. We generate a mean escape time Ti to be associated with each lattice site i using a given probability distribution. 
In this paper we use the probability distribution pir) = ^t~^~^6{t — 1). The times are fixed throughout the 
entire simulation. 

3. We distributed particles on the lattice sites with concentration c so that the probability that any particular site 
is initially occupied is c. There is at most one particle per site. 

4. The dynamics then proceeds as follows: 

(a) We generate the waiting times for each particle's next jump using an exponential distribution about a mean 
escape time Ti associated with each site i. 

(b) We choose the particle with the smallest waiting time. This particle jumps to one of its nearest neighbors 
(two in Id, six in 3ci) with equal probability. 

(c) If the destination site is empty, we update the waiting time for the arriving particle (we simply add time 
to the old waiting time a new time obtained using [3a| and we repeat the process I4bl 

(d) If the destination site is occupied, the particles annihilate (in the A + A ^ Q) or coagulate (in the case 
A-\- A ^ A). We then repeat the process by returning to I4bl In the A + A ^ A case we need to specify 
which of the two particles is the victim in the coagulation process. In the "kill" scenario the arriving 
particle (the particle that just performed the jump) annihilates the particle that was already there, and the 
waiting time of the arriving particle is updated as in[3cl that is, as if the destination site were empty. In the 
"no kill" scenario the arriving particle is annihilated, and the waiting time of the surviving particle remains 
unchanged. One could also implement a combination of these rules with some probability weighting. In 
any case, in our simulations interestingly we have observed no statistically significant difference between 
the results obtained with "kill" and with "no kill" rules. 
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